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Abstract. Three-dimensional numerical simulations with C0 5 BOLD, a new radiation hydrodynamics code, result 
in a dynamic, thermally bifurcated model of the non-magnetic chromosphere of the quiet Sun. The 3-D model 
includes the middle and low chromosphere, the photosphere, and the top of the convection zone, where acoustic 
waves are excited by convective motions. While the waves propagate upwards, they steepen into shocks, dissipate, 
and deposit their mechanical energy as heat in the chromosphere. Our numerical simulations show for the first 
time a complex 3-D structure of the chromospheric layers, formed by the interaction of shock waves. Horizontal 
temperature cross-sections of the model chromosphere exhibit a network of hot filaments and enclosed cool regions. 
The horizontal pattern evolves on short time-scales of the order of typically 20 - 25 seconds, and has spatial scales 
comparable to those of the underlying granulation. The resulting thermal bifurcation, i.e., the co-existence of 
cold and hot regions, provides temperatures high enough to produce the observed chromospheric UV emission 
and - at the same time - temperatures cold enough to allow the for mation of molecules (e.g., carbon monoxide). 
Our 3-D model corroborates the finding bv lCarlsson fc SteinI Jl994T) that the chromospheric temperature rise of 
semi-empirical models does not necessarily imply an increase in the average gas temperature but can be explained 
by the presence of substantial spatial and temporal temperature inhomogeneities. 
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• i— i , 1. Introduction A strong motivation for three-dimensional time- 
' dependent modelling arises from the need to recon- 
^ ; Three-dimensional, time-dependent radiation hydrody- ci i e a pp are ntly contradictory solar observations: Carbon 
namics simulations of solar and stellar surface convec- monoxide absorption lines imply gas temperatures as 
tion have now reached a level of sophistication which goes \ ow as ~ 3700 K in the chromosphere of the quiet 
far beyond that of idealised numerical experiments, and Sun (see iNoyes fc Halllll972t lAvres fc Testermanlll98lt 
allows a direct confrontation of such models with real [ Solanki et alJ Il994t lUitenbroek et all Il994t lllitenbroekl 
stars (e.g.. fStein fc Nordlundl Il998t lAsplund et al] |2000j |2000aHAvresll2002[ and references therein) , whereas chro- 
|Freytag et al] |200j |Ludwig et al] |200j . Extending this m0S pheric UV emission features re quire much higher tern- 
kind of simulation to include the low chromosphere, it is peratures at the sam e heights (e.g.. lAvres fc Linskvlll976T: 
possible to study - in a single model and based on first ICarlsson et al.lll997|) . 

principles -the generation of waves by the convective flow Semi-empirical models which have been con- 
as well as the wave propagation and dissipation in the struc te d based on UV and microwave observations 
higher layers. Extended simulations of this type may then , IVernazza et alJll98ll hereafter VAL; iMaltbv et alJ 
be utilised to explore the hitherto poorly understood 3-D [^jj iFontenla et all Il993i hereafter FAL) commonly 
thermal structure and dynamics of the non-magnetic chro- feature ft temperature m i nimum f Tmln « 4200 - 4400 K 
mospheric internetwork regions, and to obtain an indepen- &t ft height of z ^ 5QQ km aWe optic&1 depth un% 
dent theoretical estimate of the amount of chromospheric &nd an outwardly increas i n g temperature above. On 
heating due to acoustic waves. the other hand, m odels based on CO observations (e.g., 

IWiedemann et alJ 11994^ show a monotonic decrease of 

Send offprint requests to: wedemeyer@kis.uni-freiburg.de temperature with height. 
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These conflicting observations and the inferred repre- 
sentative models have led to a controversy about the na- 
ture of the chromosphere of the non-magnetic quiet Sun 
which is going on for many years now (see, e.g.. lKalkofenl 
l200ll) : Is the chromosphere of the average quiet Sun a 
time-dependent phenomenon with a mostly cool back- 
ground and large temperature fluctuations due to upward 
propagating shock waves? Or is it persistent and always 
hot with only small temperature fluctuations? In short: Is 
the non-magnetic solar chromosphere hot or cool? 

A large number of observations show that the chro- 
mosphere of th e quiet Sun is indeed a very dynamic phe- 
nomenon Ce.g.. ICarlsson et aljll997t iMuglach fc Schmidt] 
l200ll iKriiger et all l200lt IWunnenberg et all 12002^ 
Obviously, static one-dimensional models can only de- 
scribe selected time-averaged properties. More realistic 
modelling should therefore be time-dependent. 

Starting in the late 1960s, the pioneering work 
on 1-D time-dependent numerical models of chromo- 
spheric heating by acoustic and magneto-hydrodynamic 



waves is due to Ulmschneider and collaborator 
a long series of papers (e.g.. lUlmschneiderl 


s. In 
1971; 


Ulmschneider fc Kalkofedll977t lUlmschneider et alJ 


1978: 


Muchmore fc Ulmschneide 


rl Il985l lUlmschneider et alJ 


19871 lTJlmschneide^lT989: 


Idintz et al.lll994l). they stud- 



ied in detail the chromospheric energy balance between 
dissipation of prescribed short-period (mostly monochro- 
matic) acoustic waves and radiative emission. In their 
models, the acoustic energy flux is supplied by a piston 
acting as a lower boundary condition. Assuming that the 
generation of acoustic waves by the "turbulent" flows 
in the upper convecti on zone can be described by the 
Lighthill-Stein theory (iLighthilJ Il953 ISteirjll967l Il968t 



iMusielak et~aD 1 1994 lUlmschneider et all 119961 [19990, 

they compute dynamic chromospheric models not only 
for the Sun but also for a sample of main-sequence stars 
and giants. Based on these models, they conclude that 
the observed "b asal flux" from the chromospheres of 
late-type stars l|Schriiverl Il987t iRutten et all Il99l[) is 
fully at tributable to the dis sipation of acoustic wave 
energy I Buchholz et alJ ^98), and that the observed 
variation of chromospheric emission can be explained by 
the ad ditional heating of magn etohydrodynamic shock 
waves (lUlmschneider et al.lE00l|) . 

The detailed rad i ation hydrodynam ics simulations 
by ICarlsBon fc Steinl (|l994 Il99fil I19971 hereafter CS) 
are another prominent example of sophisticated 1-D 
time-dependent modelling. These authors successfully ex- 
plained the Ca 11 bright points as a result of propagat- 
ing shock waves. In their model, the waves are excited by a 
piston which is driven by a velocity variation derived from 
observed oscillations at the photospheric level. Instead of 
a temperature minimum and a monotonic temperature in- 
crease above, as characteristic of the VAL and FAL mod- 
els, CS find a chromosphere with a mostly cool background 
and large temperature fluctuations due to upward propa- 
gating shocks. Even more remarkable is the fact that they 
are able to reproduce the rise of the radiation temperature 



without an increase of the mean gas temperature. Basic 
reasons are the nonlinear temperature dependence of the 
Planck function in the UV and the extreme temperature 
peaks associated with the shock waves. This led CS to 
the conclusion that the chromosphere of the quiet Sun is 
not persistent but a spatially and temporally intermittent 
phenomenon which - if averaged over space and time - is 
mostly cool and not hot. 

Although the one-dimensional models of the non- 
magnetic solar chromosphere mentioned above are highly 
elaborate, including a fully time-dependent H ionisation 
and detailed NLTE radiative transfer, they suffer from the 
need for an external prescription of the wave excitation, 
and of course they cannot account for horizontal inhomo- 
geneities and the associated effect of dynamic cooling on 
the atmospheric energy balance. 

In this r egard, the three-dime nsional self-consistent 
modelling bv lSkartlien et all |2p00) can be considered as a 
major improvement. The idea of Skartlien and co-workers 
was to extend the standard radiati on hydrodynamics 
simul ations of the solar granulation l|Stein fc Nordlundl 
1998) into the chromosphere, where local thermodynamic 
equilibrium (LTE) is known to be a poor approxima- 
tion. In order t o adapt it to chromospheric conditions, 
ISkartlier] l|200fj|) upgraded the radiative transfer part of 
the Nordlund- Stein code by implementing an iterative 
method to treat coh erent isotropic s cattering in 3-D. The 
simulations enabled ISkartlien et alJ to analyse the gener- 
ation, propagation, and dissipation of acoustic waves in 
three dimensions. The main emphasis of their study was 
on the excitation of transient wave emission resulting from 
the collapse of small granules, and the dynamic response 
of the chromospheric layers to such acoustic events. 

In the present paper, we present similar time- 
dependent 3-D models which extend from the upper 
convection zone to the middle chromosphere. The ra- 
diation hydrodynamics simulations are performed with 
C0 5 BOLD, a new radiation hydro dynamics code deve l- 
oped by B. Freytag and M. Steffen l|Frevtag et alJl2002|) . 
In this exploratory simulation, we treat the radiative 
transport in LTE with grey opacities (see Sect. 12.21 and 
discussion in Sect.|SJ). This simplification allows us to work 
at a signifi cantly higher spa tial resolution (140 x 140 x 200 
cells) than ISkartlien et alJ (32 x 32 x 100 grid). We find 
that the 3-D structure of the non-magnetic chromospheric 
layers is characterised by a complex pattern of interact- 
ing shocks, forming a network of hot filaments and en- 
closed cool "bubbles". This chromospheric pattern and 
its implications are chosen as major subject of this paper 
since the topology and the dynamics of the pattern are 
likely not to be too sensitive to the LTE simplification. 
We conclude that the low chromosphere exhibits a promi- 
nent thermal bifurcation: hot and cool regions exist side 
by side. Surprisingly, this sm all-scale (non-m agnetic) net- 
work was not mentioned by ISkartlien et alJ : presumably, 
it was not noticed due to the poor (horizontal) spatial 
resolution of their numerical model. 
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In Sect. |21 we will give a short overview of the numer- 
ical details of C0 5 BOLD. The 3-D model is described in 
Sect. |3 followed by the results in Sect.0] Finally, a discus- 
sion and conclusions are presented in Sect. and Sect.[SJ 
respectively. 

2. The radiation hydrodynamics code C0 5 BOLD 

CO 5 BOLD solves the time-dependent hydrodynamic equa- 
tions coupled with the radiative transfer equation for a 
fully compressible, chemically homogeneous plasma in a 
constant gravitational field in two or three spatial di- 
mensions. Operator splitting separates Eulerian hydro- 
dynamics, 3-D tensor viscosity, and radiation transport. 
Magnetic fields are not included so far, restricting this 
version of C0 5 BOLD to internetwork regions. 

The most important properties of the c ode are de- 
scrib ed below (see also iFrevtae et all l2002t IWedemeverl 
120031) . A more detailed paper on the code itself is in prepa- 
ration (Freytag, in prep.). 

2.1. Hydrodynamics 

The relations for the conservation of mass, momentum, 
and energy are solved on a fixed Cartesian grid allow- 
ing spatially non-equidistant meshes. Directional operator 
splitting transforms the 2-D/3-D problem into 1-D sub 
steps which then can be tre ated with a fast approximate 
Riemann solver 1986). The scheme is modified to 

account for a realistic equation of state and an external 
gravity field. 

Additionally, a small amount of tensor viscosity is 
added in a separate sub step. Although the hydrodynam- 
ics scheme is stable enough to handle 1-D and most multi- 
dimensional problems, there are special multi-dimensional 
cases which require an additional tensor viscosity to en- 
sure stability. Such cases occur, e .g., near stro ng shocks 
which are aligned with the grid (Quir Our nu- 

merical scheme has proven to be very robust in handling 
shocks, which is important when modelling chromospheric 
conditions. 

2.2. Radiation transport 

The equation of radiative transfer is solved applying long 
characteristics ("rays"). A large number of rays traverse 
the computational box under different azimuthal and in- 
clination angles. Independently along each ray, the radia- 
tive transfer equation is solved with a modified Feautrier 
scheme. The radiation transport is treated in strict LTE so 
far. In this work, a grey (frequency-independent) radiation 
transport with realistic opacities is used (see Sect. 12 .3(1 . 
The applied scheme is well-suited for the lower layers 
(convection zone and photosphere), but clearly requires 
further improvements for chromospheric conditions where 
substantial deviations from LTE prevail and the UV ra- 
diative transfer is dominated by scattering. See also the 
discussion in Sect. 03 



2.3. Equation of state and opacities 

The equation of state takes into account partial ionisa- 
tion of H and He, as well as formation and dissociation 
of H2, assuming thermodynamic equilibrium. It is solved 
by interpolation in a table which is computed in advance 
for a prescribed chemical composition of hydrogen, he- 
lium, and a representative metal. The table consists of 
two-dimensional arrays as functions of density and inter- 
nal energy. 

For the model presented in this work we used 
a Rosseland mean opacity look-up table which has 
been compiled and processed based on data of 
OPA L for temperatures above 12000 K ijlglesias et alJ 
Il992l) and PHOENIX for temperatures below 12000 K 
l|Hauschildt et al.lll997l and references therein). The ta- 
ble provides the opacity as a function of temperature and 
gas pressure. 

Although a large number of atomic lines and molecular 
features are formally taken into account in the construc- 
tion of the opacity table, it is clear that the stronger lines 
are not properly represented when computing the grey 
opacity according to the Rosseland averaging procedure. 
Consequently, the stronger spectral features are essentially 
ignored in the present approach (see also Sect.EJ- 

2.4. Boundary conditions 

Located deep in the convectively unstable layers, the lower 
boundary is open, i.e., material is allowed to flow in and 
out of the computational box. The inflow of material is 
constrained to ensure a vanishing total mass flux across 
the lower boundary so that the total mass in the com- 
putational volume is preserved - aside from smaller gains 
or losses across the upper boundary. The entropy of in- 
flowing material is a prescribed parameter, and indirectly 
controls the effective temperature of a model. The vertical 
derivative of the velocity components is zero. The pressure 
in the bottom layer is kept close to plane-parallel by arti- 
ficially reducing horizontal pressure fluctuations towards 
zero with a prescribed time constant. 

At the upper transmitting boundary the vertical 
derivative of the velocity components and of the inter- 
nal energy are zero; the density is assumed to decrease 
exponentially above the top boundary. Material can flow 
into the computational box if the velocity at the boundary 
is directed downwards. The temperature of the inflowing 
material is then altered towards a temperature T top on a 
characteristic time scale of typically a few seconds. This 
simple boundary condition turns out to be stable and al- 
lows (shock) waves to leave the computational box without 
noticeable reflections. Moreover, we have chosen the loca- 
tion of the upper boundary such that it is far away from 
the regions which are of particular interest in this work. 

The lateral boundary conditions are periodic. 
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Fig. 1. Logarithmic temperature in vertical 2-D slices taken from the 3-D model at different horizontal positions 
( a) y — 1540 km, b) y — 2820 km ): Top of the convection zone, photosphere, and low/middle chromosphere with 
propagating shock waves. The solid line marks the height for optical depth unity. The dotted lines are contours for 
logT = 3.7,3.95,4.00,4.05, ..,4.20 (top to bottom). The temperature ranges from rj 16400 K to RJ 5200 K in the 
convection zone (i.e., below z = km) and decreases to RJ 3000 K in the photosphere and even down to RJ 1800 K in 
the chromosphere. 



3. The 3-D model 

The 3-D model consists of horizontally 140 grid points 
(x, y) with a constant resolution of 40 km, leading to a 
horizontal size of 5600 km which corresponds to an angle 
of rj 7". 7 in ground-based observations. The total verti- 
cal height is 3110 km, reaching from the upper convection 
zone at z = —1400 km to the middle chromosphere at a 
height of z = 1710 km. The origin of the geometric height 
scale (z = km) corresponds to the temporally and hor- 
izontally averaged Rosseland optical depth unity. In the 
following we refer to the photosphere always as the layer 
between km and 500 km in model coordinates, and to 
the chromosphere as the layer above. The 200 vertical grid 
points are non-equidistant, with a resolution of 46 km at 
the bottom which decreases with height down to a con- 
stant distance of 12 km for all layers above z = —270 km. 
The computational time step is typically 0.1 to 0.2 s. 

As an initial model, we extended an already evolved 
model which reached up to the top of the photosphere. 
The temperature and density stratification for the new 
grid cells were calculated under the assumption of hydro- 
static equilibrium. Interestingly, the further evolution of 
the model does not depend strongly on the initial con- 
dition because the chromosphere turns out to be highly 
dynamical on short time-scales. After only a few minutes 
of simulation time the initial chromosphere already formed 
the typical structures which we will discuss below. 

However, the first 170 min of the simulation sequence 
are not used for data analysis to ensure that the model has 
sufficiently relaxed. The results presented in this work are 
based on another 151 min of simulation time. 



4. Results 

4.1. Structure of the model atmosphere 

Figures ll!3l show the temperature in vertical and horizon- 
tal slices of the 3-D model which is described in Sect. 
The data for Figs. 11121 are taken from the same time step. 
Figure |21 illustrates the depth-dependence of the structure 
of the model atmosphere by means of 2-D temperature 
slices at various geometrical heights. The same figure also 
shows synthetic images of the emergent continuum inten- 
sity at A = 5500 A and A = 1600 A which were com- 
puted subsequent to the simulation for the selected time 
step. For these calculations LTE radiative transfer was as- 
sumed. We used pure continuum opacities (dominated by 
Si I b-f absorption at 1600 A) taken from the Kiel spec- 
trum synthesis package LINFOR. 

Obviously, there are striking differences between the 
horizontal patterns in the photosphere and the layers 
above. The temperature at the bottom of the photosphere 
(Fig. 0a) reveals the granulation which comes out more 
clearly in the intensity image for A = 5500 A (Fig. Elg). 
The granulation is very similar to observations in various 
aspects like shape, size distribution, and lifetime of the 
granules, indicating that in the lower p art of the model th e 
physics are realistically represented (Wedemeverl feoO:^. 
Only 250 km above, a reversed granulation pattern ap- 
pears (Fig. [2Jb): the inner parts of the granules are dark 
due to the rapid cooling of the ascending gas, and bright 
rims (note the double structure) appear at the edges of 
the granules, representing hot shocked gas being directed 
into the intergranular lanes. 

Higher up, the model chromosphere is characterised by 
a network of hot matter and small-scale hot spots on a cool 
background as can be seen in the horizontal cross-sections 
in Fig. [3d-f. The pattern is a result of interaction of prop- 
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Fig. 2. Temperature in horizontal 2-D slices at different heights in the photosphere at z — km, 250 km, and 500 km 
(a-c), and in the chromosphere at z = 750 km, fOOO km, and 1250 km (d-f). Panels g) and h) show the emergent 
continuum intensity at A = 5500 A and A = 1600 A, respectively. 
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T-tr^-iciu-i |l | Tt^hmuiII T-inpiAr+p | 

Fig. 3. Temperature in horizontal 2-D slices at z = 1000 km for a short time sequence (At — 30 s). 



agating hydrodynamic shock waves which are an ubiqui- 
tous phenomenon in the model chromosphere. The shock 
fronts are usually inclined, so a horizontal cut through the 
temperature field shows a filamentary structure. There is 
also a clear signature of oscillations with periods in the 
3-min range (see Fig. UJ. Shock waves are present at all 
time steps, mostly several at the same time (Tigs. ITl3*|l . 
The waves propagate in the vertical as well as in the hori- 
zontal direction and interfere with each other, compressing 
and heating the gas in the filaments (see Sect. I4.2|l . 

As a consequence of the correlation between convective 
motions and the excitation of acoustic waves, the spatial 
scale of the pattern is comparable to that of the underlying 
granulation. 

The network-like pattern appears more subtle in the 
UV continuum intensity at a wavelength of A = 1600 A 
(Fig.[2h). Rather, a small area of enhanced emission 
stands out of an otherwise dark background. This is 
caused by the highly non-linear temperature response of 
the Planck function in the UV. Hence, the hot gas, which 
is connected to the propagating shock waves, contributes 
by far more to the emergent UV continuum intensity than 
the cool regions. Note that for more realistic results, scat- 
tering and line blocking must be taken into account. 

Due to the ongoing propagation of the waves the pat- 
tern changes continuously (see Fig.^J) on time-scales which 
are much shorter than derived for the granulation. We 
calculated autocorrelation times for sequences of horizon- 
tal temperature slices and determined height-dependent 
pattern evolution time scales as the time lags for which 
the autocorrelation decreased to a value of 1/e. At chro- 
mospheric heights the characteristic time scales are as 
short as 20 — 25 s whereas the same analysis produces 
time scales of ^ 120 s at the bottom of the photosphere 
(z = 0). Using the emergent grey intensity, which ren- 
ders the low photosphere, instead of the gas temperature 
leads to ~ 200 s. The difference between temperature 
and intensity result can be understood if one considers 
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Fig. 4. Variation of relative velocity power with height. At 
each height the power spectrum of the horizontal average 
of the vertical velocity was integrated over the frequency 
intervals which are specified on the left. 



that structures also move up and down, for instance, due 
to oscillations. Consequently, the pattern at a fixed geo- 
metrical height changes more quickly than visible in the 
corresponding intensity. Furthermore, spatial smearing of 
the pattern, i.e., reducing the image resolution to values 
caused by observational seeing conditions, produces longer 
time scales. This should be kept in mind when comparing 
the theoretical results with empirical data. 

4.2. Waves, oscillations, and shocks 

Acoustic waves are excited by various processes con- 
centrated in the uppermost layers of the solar convec- 
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k p n\) k j I m| • f n$ 

Fig. 5. Formation and propagation of shock fronts: Each column shows a time sequence of vertical slices (temperature) 
taken from different positions and times of the 3-D model, a) (left column) arch- like/spherical wave, b) (middle 
column) plane wave, c) (right column) waves excited by merging downdrafts. Note the different time steps which are 
quoted in the lower left corners. The temperature is colour-coded for the range T = 2000 K to T = 7500 K. Additional 
contour lines are present for T = 5000 K (dotted), and T = 7500 K, 9000 K, and 10000 K (all solid). 



tion zone. Excitation processes have be en investigated 
by me a ns of hydrodynamical mod elling bvlSkartlien et aT 
2QQ3), iNordlund fc Stehl JUS), and IStein fc NordW 



2001). Skartlien et al. study the collapse of small gran- 



ules which leads to transient wave emission. Nordlund & 
Stein focus on the interaction of convection with resonant 
oscillatory modes to derive an estimate of the power input 



into the solar 5 min oscillations. Like the afore mentioned 
authors, we observe in our model the excitation of both 
propagating and standing acoustic waves. The standing- 
waves are the model analogs to the solar 5 min oscilla- 
tions. Together with the propagating waves they generate 
a complex interference pattern in the photospheric and 
chromospheric layers, where shocks are frequently formed. 
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Fig. 0] illustrates the distribution of power among ra- 
dial oscillations as a function of height. Fourier spectra 
were calculated for a 151 min long time sequences of 
the horizontally averaged vertical velocity component at 
each height independently, and integrated over frequency 
bands roughly centred around periods of 5 min and 3 min. 
Figure 01 shows that the dominant contribution to the ve- 
locity power shifts from the 5-min band to the 3-min band 
at around z ~ 1200 km. We find no significant power in 
the low frequency band (periods larger than ~ 420 s), 
while the high frequency band (periods below ~ 140 s) 
contributes some power in the higher layers. 

The absolute energy of the oscillatory motions (not 
shown) decreases in all bands with increasing height. The 
largest energies are found in the deepest layers, indicat- 
ing that the excitation of the oscillations takes place in 
the convection zone for all frequencies. The 5-min band 
lies below the acoustic cut-off frequency (~ 5.5 mHz) ren- 
dering these waves evanescent while in the 3-min band 
some frequencies allow propagating waves. This implies a 
stronger damping in the 5-min band, and explains why 
the "3-min" oscillations dominate in the chromosphere: 
the decline of energy with height is more pronounced in 
the 5-min band than in the 3-min band. A localised non- 
linear process converting oscillatory energy in the 5-min 
band into energy in the 3-min band is not readily appar- 
ent. 

Examples of propagating waves and shock formation 
are shown in Fig. [5] The example of the left-most col- 
umn (Fig. 0a) is displayed more quantitatively in Fig. 
for further discussion below. It shows the case of a rather 
localised shock which was triggered by pressure distur- 
bances emerging from the downflow region visible in the 
deeper layers. The formation of a spherically shaped shock 
is a frequent pattern. The spherical shock front appears as 
an upward travelling arch-like feature in our 2-D cuts. The 
middle column in Fig.0shows an example of a front which 
is horizontally more extended. In movies such events ap- 
pear often as if the front detaches over a broader area from 
the photospheric granulation pattern. It can extend over 
more than one granule and tends to preserve the shape of 
the granular pattern for some time. In the simulation, pref- 
erentially resonant modes of long horizontal wavelength 
are excited. They provide the horizontally coherent os- 
cillations which are necessary to produce these extended 
horizontal wave fronts. The right-most column (Fig. [5Jc) 
shows the formation of shocks above merging downdrafts, 
i.e., downflows in the intergranular lanes. This kind of 
event corresponds to the collapse of sma ll granules and 
has al ready been investigated in detail bv lSkartlien et alJ 
From the vertical slices in Fig. 0c it can be seen 
how two downdrafts are advected horizontally and even- 
tually merge, producing a stronger and more extended 
downdraft. During the process upward propagating waves 
are excited which may transform into shocks in higher 
layers. Moreover, a strong downdraft is often accompa- 
nied by shocks of a different nature. They come about by 
fast horizontal flows towards the downdraft. Shocks form 
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Fig. 6. Vertical profiles of the flow field shown in Fig. 0a 
at different times along the horizontal position x = km. 
We plot the temperature (solid) , the vertical velocity com- 
ponent (triple-dot-dashed), and the logarithmic pressure 
(dashed) on a linear scale. The data range to 1 corre- 
sponds to K to 7000 K in temperature, —1.59 to 4.50 
(cgs units) in the logarithmic pressure, and kms -1 to 
15 kms -1 in the vertical velocity component. The vertical 
grid is shown at the top of the figure. 



where the flow is turned into the downdraft. In Fig. 0c 
they are visible as roughly vertical features attached to 
the edges of a downdraft. These shocks interact with the 
shocks associated with the wave field (see frames at 60 s to 
120 s). Note that Fig.0shows particularly clean examples 
of the types of shock events encountered in the simulation. 
Usually, the pattern of shocks is very entangled, and often 
all features discussed before are present at the same time. 
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Fig. 7. Temperature histograms for the 3-D model (151 min of simulation time). For each height step (and all time 
steps) the temperature values of all grid cells within the corresponding horizontal plane are sorted into temperature 
bins of AT — 100 K. The vertical axis denotes the relative abundance of grid cells within a temperature bin with 
respect to all cells at that height, a) Height-dependent histogram surface, b-c) Histograms at fixed heights in the 
photosphere and d-e) in the chromosphere. The dotted lines represent the corresponding mean temperature. 



The wave depicted in Fig. (see also Fig. 0a) is an 
extreme example as a positive vertical velocity of v z « 
11 kms -1 is reached in the chromosphere. Most velocities 
are smaller. We find approximate upper limits for 95 % 
of all upward directed vertical velocities, depending on 
height: ra 4.9 kms^ 1 at z = 800 km and ~ 7.0 kms -1 
at z = 1000 km. In contrast to one-dimensional simula- 
tions, the waves in our 3-D model do not only propagate 
in the vertical direction but also horizontally. At a height 
of z = 1000 km we find that 95 % of all grid cells exhibit 
horizontal velocities of less than w 12 kms" 1 and 50% 
have values of w 5 kms -1 and below. 

An important point is illustrated in Fig. shocks are 
preferentially formed in low-density material which is flow- 
ing down from above at high velocities. The material has 
been pushed upwards by a precursory wave and now falls 
back again. The shock front is travelling upstream into the 
down-flowing material. In extreme cases the downflowing 
material is close to free-fall conditions, and flow veloci- 
ties exceed the lo c al sou nd speed. The 1-D simulations by 
ICarlsson fc Steinl lll997h exhibit a similar shock structure 
(see their Fig. 14). Judging from the same figure, Carlsson 
& Stein find typically at most one well developed shock 



in the photospheric and chromospheric layers at any given 
instant in time. Looking at one particular vertical column 
in our 3-D model we make a similar observation, finding 
typically one, sometimes two fronts. While in their piston- 
driven model Carlsson & Stein derive the wave excitation 
semi-empirically from observed time sequences of photo- 
spheric oscillations, the shock frequency in our case is a 
natural outcome of the simulation. The spatial shock fre- 
quency translates into a temporal recurrence of shocks on 
a time scale of ~2-3 min (see also Fig. 110(1 . 

4.3. Thermal bifurcation 

Although the chromospheric pattern evolves on very short 
time scales (see Sect. l4.l|) . the general picture remains the 
same in time, i.e., the chromosphere appears as a network 
of hot matter with intermittent cool regions. This ther- 
mal bifurcation can be quantified via a height-dependent 
temperature histogram. For each horizontal slice in the 
model (constant height z for each slice) a histogram of 
the temperature values is calculated for temperature bins 
of AT = 100 K. The result is shown in Fig. CI In the pho- 
tosphere, the temperature is distributed close to a mean 




Fig. 8. Temperature stratifications of different models on a geometric height scale (a) and on an optical depth scale 
(b): Horizon tally and tempo r ally av eraged grey e missivity temperature a nd mean gas temperature for the 3-D model, 
model C bv IVernazza et all i|l98lj) . mo del A bv iFontenla et aD £l993), mean g as temperatu r e and semi-empirical 
stratification of the dynamical model bv lCarlsson fe StemT lIlQQffl . and COOLC bv lAvres et aD (Il98fj) . 



value with only moderate deviations, whereas in the chro- 
mosphere, the distribution splits up into low and high tem- 
peratures. Again, this indicates the co-existence of a cool 
background and hot shocked material. 

To facilitate a ro ugh comparison with multi- 
component models (e.g.. lAvres et al~lll986l lAvrettlH995t 
lAvres fc Rabmlll99ri: lAvresl 120021) " we give approximate 
values for a hot and a cool component of our model 
chromosphere (intermediate values are neglected): Above 
z = 800 km the hot temperature ridge in Fig. [7\ peaks 
at Thot = 5500 — 5900 K, whereas the cool temperature 
peak decreases with height from T coo i = 2600 K at 
z = 800 km to » 2000 K for the upper layers of the 
model chromosphere. Thus, the hot component is com- 
parable to the temper atures in the s e mi-em pirical models 
C by VAL and C bv iMaltbv et~al] l|l986h in the height 
range 800 - 1000 km and 900 - 1100 km for the model 
A by FAL, respectively . The cool c ompo nent is much 
col der than COOLC bvlAvres et al] lll98fih and COOL0 
bvlAvres fc Rabml ljl99fih . It is much more like COOL1 
llAvresll2002r around z = 800 km which, however, is only 
valid if there is a dominating warm component. 



4.4. Temperature stratification 

In this section we discuss the consequences of the ther- 
mal bifurcation for the average temperature stratification. 
The horizontally and temporally averaged gas tempera- 
ture for the sequence of 151 min simulation time from our 



model (thin solid line in Fig. |B]a) decreases with height 
until it reaches values between 3800 K and 3700 K above 
z ~ 730 km, i.e., in the chromosphere. It does not show a 
notable temperature minimum nor a significant tempera- 
ture increase in the chromosphere like it is the case in the 
semi-empirical models by VAL and FAL (see Fig.[H]a). 
This is qualitatively similar to t he mean gas temp e rature 
profile in the 1-D simulation bv ICarlsson fc Steinl l(l995h 
which also does not show a temperature increase (see 
Fig. |5]a). However, we obtain chromospheric gas temper- 
atures which are much lower than in the simulations by 
CS. In fact, our mean chromospheric gas temperature lies 
about 1000 K below the (grey) radiative equilibrium tem- 
perature of 4680 K. The mean temperature stratification 
is rou ghly comparable to model COOLC bv lAvres et alJ 
( 1986|), which was constructed as the cool constituent in a 
multi-component model (see Fig. [H]a, where we converted 
the original column mass density scale into a geometrical 
height scale on the basis of model C by VAL). 

The semi-empirical models are based on spatially and 
temporally averaged intensities and thus refer to a static 
and homogeneous chromosphere. We note that the mean 
gas temperature from our model matches almost perfectly 
the semi-empirical models up to a height of z « 500 km. 
Above that height, the thermal bifurcation becomes in- 
creasingly significant, i.e., the temperature fluctuations 
become large (see Figs. and EJ). Clearly, the assumption 
of spatial and temporal homogeneity is not valid in the 
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chromosphere, and any one-dimensional static description 
must fail. 

CS pointed out that the chromospheric temperature 
rise in the semi-empirical models is only an artifact caused 
by the "temporal averaging of the highly nonlinear UV 
Planck function" . Furthermore, CS confirmed this by cal- 
culating a temperature distribution for their dynamical 
model in a similar way as VAL. They adjusted a steady- 
state temperature stratification to reproduce the time- 
averaged continuum and line intensities as a function of 
wavelength which are a result of their dynamic simulation. 
The semi-empirical model derived in this way by CS is a 
much better fit to the models VAL and FAL (see Fig. |BJa). 

Since no wavelength-dependent intensities are avail- 
able for our simulation (except for a few images similar to 
those shown in Fig. 0g-h), we calculated a qualitatively 
similar quantity, namely an "average grey emissivity tem- 
perature", by averaging the grey emissivity KpT 4 , where 
n is the opacity and p the density. The corresponding emis- 
sivity temperature T cm is then evaluated as: 



' (KpT 4 ) x . y 
, ( K p)x,y 



1/4 



(1) 



The brackets (-)x,y, { •}* indicate horizontal and tem- 
poral averaging, respectively. The resulting average 
temperature profile, calculated on a geometrical scale 
(thick solid line in Fig. |SJa) is indeed similar to model C 
by VAL and model A by FAL. It exhibits a temperature 
minimum at approximately the same height; the tempera- 
ture values reached in the middle chromosphere are com- 
parable. This qualitative match is better than expected 
from such a crude approximation. Thus, like CS we are 
able to produce an emissivity temperature stratification 
qualitatively similar to the semi-empirical models, with- 
out a significant increase in the mean gas temperature. 

The averages presented so far are calculated on a 
geometrical height scale. In contrast, the average grey 
emissivity temperature and the simple arithmetic average 
shown in Fig. |S|b are calculated on the Rosseland optical 
depth scale which already incorporates the distribution 
of opacity and density. Hence, the emissivity temperature 
is given by T 4 averaged over surfaces of constant optical 
depth. 

We note that the mean chromospheric gas temperature 
obtained from averaging on the optical depth scale are 
systematically higher (but still below the radiative equi- 
librium value) than those on the geometrical height scale; 
the minimum values differ by more than 500 K. That is 
caused by the fact that fluctuations appear much smaller 
on sur faces of equal optical depth (see e.g., lUitenbroekl 
l2000b|) . In a wave front the optical depth increases signif- 
icantly. Thus, averaging on an optical depth scale is done 
on surfaces which are not plane but shaped by the spatial 
inhomogeneities while averaging on a geometrical height 
scale is done on strictly plane surfaces which cut through 
the inhomogeneities. Consequently, the temperature dis- 
tribution on a surface for a particular optical depth dif- 
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Fig. 9. Horizontally and temporally averaged temperature 
fluctuation: Absolute deviations dT lms (solid, left axis) 
and relative deviations dT rms /To (dashed, right axis) on 
the geometrical height scale (a) and on the optical depth 
scale (b). 

fers from the one for a corresponding geometrical height, 
leading to different horizontal averages and thus different 
temperature stratifications. 

4.5. RMS-temperature fluctuations 

Here, we quantify the rms-temperature fluctuations which 
are another measure characterising the thermal structure. 
They arc defined by 



To To 



(2) 



where T = (T) x,y,t is the temporally and horizontally av- 
eraged temperature stratification. The quantity dT rms /To 
has been calculated on a geometrical and on an optical 
depth scale for the same model sequence as in Sect. 14.41 
(Fig. EJ). It is strongly height-dependent as can also be 
seen directly from the horizontal slices in Fig. [3 for differ- 
ent heights and from the temporal temperature variation 
in Fig. EH Obviously, the lower layers of the solar atmo- 
sphere in our model are relatively homogenous with only 
small temperature fluctuations, in contrast to the inho- 
mogeneous chromosphere. 

Like for the temperature stratification (Sect. ^3} there 
is a difference between the geometrical height scale and 
the optical depth scale. Again the temperature deviations 
are generally much smaller on a surface of a particular 
optical depth than for a corresponding geometrical height 
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Fig. 11. Height-dependent duration of cool episodes for 
a threshold temperature of Tthrcs = 3700 K in the 
3-D model, a) Absolute values for the average (solid) 
and the average plus/minus standard deviation (dotted); 
b) Ratio of integrated cool time to total time. 



(see, e.g. . IUitenbroekl2000rj|) . In both cases the average lies 
below dT lms /T rj 0.42. For particular vertical positions 
and time steps, maximum values of ~ 1.0 can be reached. 

A compara ble quantity ST/T has been used by 
iKalkofenl ll200lft to distinguish between the two opposing 
cases of a hot chromosphere with small temperature fluc- 
tuations (ST/T « 0.1) and a cool one with large fluctua- 
tions (ST/T « 10). Our model lies in between these cases. 
As mentioned earlier, the inclusion of time-dependent ion- 
isation likely leads to higher temperature peaks and ac- 
cordingly to larger temperature deviations. 

4.6. Cool regions 

As a consequence of the propagating shock waves, the tem- 
perature at a fixed position in the model chromosphere 
varies by several 1000 K with time, featuring sharp tem- 
perature peaks on top of a cool background (Fig. 1101c). 
Observations of the infrared CO fundamental vibration- 
rotation lines imp ly temperatures as low as 3700 K (e.g, 
IUitenbroekll2000a|) which also represent a lower limit for 
the average temperature stratification of the 3-D model 
(see Sect. 14.40 . We adopt this temperature as a thresh- 
old value and determine how long the temperature at 
a fixed position in the model stays below this value. In 
the following we will refer to these time intervals with 
T < Ttiucs = 3700 K as cool episodes. In Fig.HHIc such 
episodes are illustrated. The duration of a cool episode is 



influenced by the local background temperature and the 
temperature fluctuations due to the propagating waves. 
Therefore, it depends on height. However, the average 
duration stays more or less constant throughout a wide 
height range in the chromosphere. For the 3-D model, we 
determined the average duration of the cool episodes in the 
chromosphere to be 70 — 100 s fFig. II Hal . In some cases 
the cool episodes are much longer, up to several hundred 
seconds. 

With regard to a more global view of the chromosphere 
not only the duration of single cool episodes but also the 
sum of all durations is interesting. The temperature at a 
fixed position in the chromosphere of the 3-D model stays 
roughly half of the time below T thros = 3700 K (Fig.EUb). 
In the lower photosphere cool episodes are rare and thus 
negligible with regard to the total time. 

The spatial scales of the cool regions might also be 
interesting for the interpretation of observations. The av- 
erage radius of a cool region is hard to determine because 
the regions are often not closed structures like a cloud but 
are connected to other cool regions in a complicated way. 
As can be seen from Fig. the spatial scales are on aver- 
age comparable to the granulation, except for some rare 
cases with larger cool areas. The fraction of the integrated 
cool area at a particular height shows only relatively small 
temporal fluctuations (Fig. ^Ja). Thus, the model chro- 
mosphere is never completely cold and never completely 
hot. There are always cool regions next to a hot compo- 
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Fig. 12. Ratio of cool area to total area for a thresh- 
old temperature of T t h TO s = 3700 K in the 3-D model, 
a) Variation with time (solid) and time-averages (dotted) 
at different heights; b) Height-dependent time-average 
(solid), ± standard deviation (dotted), and maximum and 
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nent. The height-dependent time- average of the cool area 
fraction (see Fig.^lb) is equal to the average ratio of cool 
time to total time (Fig. ^Jb) because both represent the 
horizontally and temporally averaged number of grid cells 
with temperatures below the threshold value. On average 
50 to 60 % of the whole time and of the whole area in a 
horizontal slice of the model chromosphere has a temper- 
ature below 3700 K. Consequently, this cool component is 
not just a minor constituent in our 3-D model. 

4.7. Carbon monoxide 

It is not obvious how the variable hydrodynamic condi- 
tions affect the formation, dissociation, and spatial distri- 
bution of CO molecules in the outer solar layers. Here we 
present the results of a simple time-dependent calculation 
of the CO concentration, demonstrating that the predicted 
height distribution of CO can be very different in a static 
and in a dynamic solar atmosphere, because the reaction 
rates are highly non-linear functions of temperature. 

For simplicity, we assume that CO is formed by direct 
radiative association, C + O — > CO + hv : and is destroyed 
by collisional dissociation, CO + H — > C + O + H. In 
this case, the temporal evolution of the CO concentration 
[CO] is governed by the differential equation 



where [CO] = nco/( n c + n co)] a value of 1 means 
that all carbon is bou nd in CO molecules. According to 
lAvres fc Rabinl 1 19961) . the constants fci and k 2 depend on 
the number density of neutral hydrogen, nu, and on the 
temperature T as 



ki = 2.510~ 5 n 15 f' - 6 
and 



fcl 1 + 40T 



22.2 



(4) 



(5) 



with the notations nis = ?i,h/(10 15 cm 3 ) and T = 
77(5000 K). 

In a static environment, the equilibrium CO concen- 
tration, 



[CO] cq = ^ = (l + 40f 22 - 2 ) 1 



is approached with a characteristic time scale 
410 4 f-° 6 



chem 
CO 



— k 2 — 



n X5 1 + 40T 22 - 2 



(6) 



(7) 



dt 



([CO]) = Ai-fc [CO], 



(3) 



The characteristic time scale i^o™ according to Eq. (7J 
and the equilibrium CO concentration [CO] oq according 
to Eq. © are plotted as a function of height for the mean 
temperature and density structure of our 3-D simulation 
in Fig. E| (thin solid lines), ^q" 1 varies by orders of mag- 
nitude, from w 0.1 s at z = km (t w 1) to ^ 10 6 s at 
z = 1000 km. This variation is partly due to the temper- 
ature dependence of ico" 1 ' but mainly due to the den- 
sity factor. In the lower chromosphere (z > 600 km), 
[CO] J; 0.9, implying that almost all carbon in the chro- 
mosphere is bound in CO. 

The situation is quite different in a dynamic atmo- 
sphere. We have investigated the time-dependent case 
quantitatively, using 196 representative vertical columns 
of the 3-D model described before. The simulations 
then provide the temperature and density variations at 
each point for the time interval of 151 min. These pre- 
scribed fluctuations translate into time-dependent coef- 
ficients fci and k 2 according to Eqns. iQJ) and ©. We 
have solved Eq. © with these time-dependent coefficients 
for each point in the selected columns, using a standard 
fourth-order Runge-Kutta method. If necessary, the time- 
sequences from the simulations are repeated until a (dy- 
namic) equilibrium is obtained. 

The resulting horizontally and temporally averaged 
CO concentration is shown in Fig. Elb (thick solid line). 
It differs dramatically from the distribution found in the 
static mean atmosphere (thin solid line), except for the 
deep photosphere (z < 200 km). In the dynamical at- 
mosphere, CO is present in the photosphere and in the 
low chromosphere with a maximum concentration of only 
<[CO]) t « 0.10 at z « 340 km. Very little CO is found in 
the layers above z « 700 km. 

This finding is i n lin e with the calculations by 
lAsensio Ramos et all (|2003l) which are based on the 1-D 
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Fig. 13. Formation and destruction of carbon monoxide, a) Chemical time scale i^Q 111 for the averaged 3-D stratification 
(solid) and for the time-dependent calc ulation t c ^ m ' dyn (dashed). The black dot indicates a time scale of 10 hours 
at the classical temperature minimum llAvres fc Rabinl Il996l) . b) Equilibr ium CO concentration for the averaged 
3-D stratification (thin solid), and average CO concentration of the time-dependent calculation (thick solid) together 
with the corresponding maximum and minimum values (dashed). [CO] = 1 means all carbon is bound in CO molecules. 



numerical simulations by CS. They, too, state that no sig- 
nificant CO concentration should be present at heights 
greater than « 700 km. Furthermore, the difference in 
the CO concentration between the static and the dy- 
namic approa ch becomes larger with in creasing height 
(see Fig. 2 in lAsensio Ramos et al.l Eo03') which is qual- 
itatively si milar to our results (s ee Fig.Elb). Hence, we 
agree with lAsensio Ramos et alJ that detailed nonequilib- 
rium CO chemistry must be taken into account. 

The reason for the difference between the static and 
the dynamic model presented in this work, however, is re- 
lated to the fact that CO is rather efficiently destroyed 
during the passage of high-temperature regions (shock 
fronts), where the chemical reaction time scales are short. 
In the subsequent cool phases (see also Sect. 14.6)1 . reac- 
tion time scales are much longer, so the CO concentration 
builds up rather slowly and reaches only moderate levels 
before the next high-temperature event occurs. The low 
CO concentration in the upper atmosphere is thus a con- 
sequence of the onset of shock formation and the related 
higher temperature peaks (see Fig. [7]a). 

If chemical time scales are short compared to the hy- 
drodynamical time scales, i^o" 1 <C £hd, then chemical 
equilibrium is reached instantaneously, and [CO]^ n « 
(fci/fe)t- Hence, we find that the CO concentration in 
the lower photosphere (z ^ 200 km) is reasonably well 
represented by this approximation. In these layers, the 



spatial CO distribution is tightly correlated with the lo- 
cal temperature: the coolest regions have the highest CO 
concentration. 

For the higher layers (z <; 200 km), which are of more 
interest in this investigation, we find t^o* 1 > inD- The 
CO concentration is then well approximated by [CO]^° « 
(ki)t/ (lt2)t- Since &2 is a highly non-linear function of T, 
the high-temperature events vastly dominate the time av- 
erage and the resulting CO concentration is much smaller 
than implied by the mean temperature. 

We can also conclude that the dynamical equilibrium 
CO concentration is attained on a characteristic time 
scale £pQ ln,dyn — (fc 2 )^ 1 , which is also shown in Fig. 1131a 
(dashed). The correlation between [CO] and T is expected 
to be poor in the higher layers: the highest concentrations 
build up in places with the longest history of relatively 
undisturbed conditions ([CO] w 0.4), and almost no CO 
is found just behind strong shock fronts. 

The results described above can only be a first esti- 
mate of the height profile of [CO] under time dependent 
conditions, because the underlying calculations still have 
severe limitations. More secure conclusions about the CO 
distribution in the upper solar atmosphere have to wait 
for more detailed future simulations taking into account 
(i) the transport of CO molecules with the flow, (ii) a more 
complete chemical reaction network incl uding multi-step 
reactions affecting the CO balance (see lAvres fc Rabinl 
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Il996t IXsensio Ramos et al.ll2003|) . and (iii) the back re- 
action of the CO concentration on the radiative cooling 
rate. 

5. Discussion 

Quantitatively, the results presented in the pre- 
ceding sections must be considered as preliminary, 
since the physics of C0 5 BOLD are not yet prop- 
erly adapted to chromospheric conditions. In partic- 
ular, the assumption of LT E is a poor approxi- 
mation in the chromosphere jCarlsson fc Steinl [2002: 
iRammacher fc Ulmschneiderll2003l) . A realistic treatment 
should account for deviations from LTE and also requires 
the time-dependent computation of the ionisation of hy- 
drogen and other important species. Nevertheless, we be- 
lieve that some of the basic features seen in our model are 
insensitive to the detailed treatment of thermodynamics 
and radiative transfer. 

The 3-D topology of the small-scale chromospheric net- 
work we discovered in our simulation, and its spatial and 
temporal scales are expected to be a robust feature. This is 
confirmed by test calculations with different values of the 
tensor viscosity, a different grey opacity table, and even 
with a frequency-dependent (multi-group) radiative trans- 
fer scheme using five opacity bins: the dynamical proper- 
ties of these models (like the height-dependent amplitude 
of the velocity fluctuations) turn out to be quite insuscep- 
tible to changes of the analysed numerical parameters. We 
attribute this to the fact that the dynamics of the model 
chromosphere are governed by the lower layers where the 
excitation of acoustic waves takes place and that the nu- 
merical modelling of these layers, i.e., the photosphere and 
the top of the convection zone, is quite realistic. In this 
context, it is reassuring to find prominent chromospheric 
oscillations in the 3-min range whose properties are largely 
independent of the numerical details of the simulation. 
The qualitative similarity to observations indicates that 
the dynamics are indeed modelled reasonably well. 

The horizontal structure of our model chromosphere, 
i.e., its topology, is reminiscent of observed patterns 
like the chromosph eric "background pattern" found by 
feriiger et alJb oOl) and the structure of Can H obser- 
vations l)Sutterlml l2003) . The latter, for instance, exhibits 
spatial scales which are comparable to the granulation 
and thus to the scales which are found in our numerical 
simulation. However, the observed patterns originate pre- 
dominantly from lower layers and should therefore not be 
confused with the patterning of the model chromosphere. 
The differences might be revealed by determining the time 
scales on which the different patterns evolve. This issue 
needs to be investigated more properly in the future. 

In contrast to the spatial scales and the topology of the 
atmospheric patterns, the amplitude of the temperature 
fluctuations in the model chromosphere is more suscepti- 
ble to the treatment of radiative transfer. Indeed, the tem- 
perature fluctuations are significantly smaller in the afore- 
mentioned test calculation using a frequency-dependent 



radiative transfer scheme. However, we recall that, up to 
the mid-chromosphcric layers (z — 1000 km), the peak 
shock temperatures i n our grey simulat i on are very simi- 
lar to those found by ICarlsson Steinl JTflol rmfll and 
by ISkartlienl l)l998fl . More precisely, in the lower and mid 
chromospheric regions CS find peak temperatures which 
lie only somewhat (^1000 K) above our values. 

The shock peak temperatures are of importance since 
many spectral features are biased towards high temper- 
atures. Theoretically, the peak temperatures depend on 
the shock strength and are given by the Rankine-Hugoniot 
jump conditions. In the absence of radiation, a conserva- 
tive numerical scheme guarantees that the jump condi- 
tions arc fulfilled, i.e., the post-shock temperature is inde- 
pendent of the spatial resolution. If radiation is important, 
however, the peak post-shock temperature is reduced rel- 
ative to the theoretical value by an amount that depends 
on the spatial resolution of the numerical grid. This is 
because the shock heating is stretched out over a finite 
time interval, given by the time a volume element needs to 
cross the shock front which is smeared out over a number 
of grid points. In this case, radiative cooling can reduce 
the attainable peak temperature if the radiative cooling 
time is comparable or smaller than the time scale of shock 
heating. Furthermore, the overall energy dissipation in the 
shock is altered due to a change of the effective adiabatic 
exponent of the gas. 

In our model based on grey radiative transfer all chro- 
mospheric layers are optically thin. Here, radiative cooling 
times are independent of the flow geometry and mainly 
dependent on temperature. The cooling time at a tem- 
perature of 7000 K — about the highest temperature we 
observe in the simulation — amounts to ^200 s, and is 
increasing rapidly for lower temperatures. The dissipation 
time scale in the shocks is in the order of a few seconds. 
This means that the thermal structure of our shocks is 
hardly affected by radiation and primarily given by the 
shock strength. Only in the most extreme cases we ex- 
pect some limiting influence of radiative cooling on the 
post-shock temperature. Similarly, the thermal structure 
of the post-shock regions is mainly controlled by cooling 
via adiabatic expansion. 

As mentioned above, the differences in the peak tem- 
peratures of our model and the simulation by CS become 
larger in the higher layers. First, CS employ an adap- 
tive grid in their simulation with a grid spacing of typ- 
ically 200 m near shocks which is thus much finer than 
our fixed (vertical) spacing of 12 km. Furthermore, it 
appears plausible that our radiative transfer - based on 
Rosseland opacities including lines as true absorption 
produces shorter radiative cooling times compared to CS. 
The higher resolution and longer radiative cooling times 
in the model of CS lead us to expect that their shock peak 
temperatures are also largely unaffected by radiation. 

Two effects can explain the somewhat higher shock 
temperatures of CS. First, the shock strength in the CS 
model might simply be higher than in our case. This could 
be related to the semi-empirical piston velocity CS feed 
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in at the bottom of their model, or their 1-D geometry- 
forcing shocks to remain plane-parallel. As we have seen 
above, extended horizontal shocks are more the exception 
than the rule in our 3-D simulation; most shock fronts 
weaken as they propagate radially away from their source. 
Another effect is related to our assumption that thermody- 
namic equilibrium conditions prevail in the chromosphere. 
In a recent paper, ICarlsson fc Steinl l|2002() demonstrated 
that this is a poor approximation. Ionisation equilibria 
cannot follow the rapid thermodynamic changes intro- 
duced by the flow. One consequence is that the energy 
which is dissipated in shocks cannot go into ionisation but 
has to go into a temperature increase of the post-shock 
gas. Moreover, accounting for finite recombination time 
scales instead of assuming ionisation equilibrium could re- 
duce the ability of the post-shock gas to cool, thus lead- 
ing to even higher temperatures. Since CS account for the 
thermodynamic non-equilibrium effects, their shock tem- 
peratures should be higher. 

We cannot decide on the basis of the available in- 
formation which is the reason for the differences in the 
peak temperatures of the shocks. However, we conclude 
that the differences up to the mid-chromospheric layers 
(z = 1000 km) are modest: our peak temperatures are 
w 7000 K compared to w 8000 K in the CS model. We 
further n ote that the chr omospheric peak temperatures 
found by ISkartlienl lll998h (see his Fig. 9) are SS 7000 K 
in their case of frequency-dependent radiative transfer ac- 
counting for line scattering, which is surprisingly close to 
our result obtained with our grey LTE radiative transfer. 

This supports our conclusion that our grey radiative 
transfer employed in this work is more realistic than the 
frequency-dependent method available for C0 5 BOLD. The 
latter method strongly overestimates the (LTE) cooling in 
the strong spectral lines (treated as true absorption) and 
thus wrongly reduces the maximum attainable tempera- 
tures. However, note that also the grey radiative transfer 
is not appropriate for chromospheric conditions. Rather, a 
detailed frequency-dependent non-LTE radiative transfer 
is necessary. Furthermore, for a quantitative comparison 
with the observations it would be necessary to perform 
three-dimensional spectrum synthesis, which is planned 
for the future. 

In contrast to the peak temperatures, the mean chro- 
mospheric temperature (and also the minimum tempera- 
tures) in our simulation are significantly lower than those 
foun d by CS ( see Fi g.lBl), and also somewhat cooler than in 
the lSkartlienl l|l998h model. Obviously, the mean tempera- 
ture structure is more strongly influenced by the treatment 
of radiative transfer than it is the case for the peak tem- 
peratures. It should therefore be considered as uncertain. 
Somewhat surprisingly, however, we note that our grey 
and our frequency-dependent simulations produce almost 
identical mean chromospheric temperature structures. 

Nevertheless, the differences in the average tempera- 
ture stratification between the one-dimensional simulation 
by CS and the presented 3-D model can be understood 
if one watches the velocity field of a region which just 



has been traversed by a strong shock wave. The flows are 
mostly directed outwards away from the centre of such a 
region. We interpret this as fast and thus adiabatic ex- 
pansion of the traversed region. This "dynamic cooling" 
is obviously more efficient in 3-D than 1-D simply due 
to the additional spatial dimensions. This effect thus pro- 
duces lower average and minimum temperatures in our 
simulations compared to those of CS. 

Furthermore, we point out that the thermal bifurca- 
tion in our 3-D model (see Sect. 14.3(1 is not due to the 
action of carbon monoxide as a cooling agent. Rather, it 
is caused by the acoustic wave field and the resulting dy- 
namic cooling of adiabatically expanding regions as dis- 
cussed above. Carbon monoxide is only taken into account 
in the grey opacity tables so far, and so its real influ- 
ence is underestimated. A similar simulation with a dif- 
ferent grey opacity t able without m olecular contributions 
(based on ATLAS6, iKuruc d Il97(i leads to very similar 
results. We thus conclude that CO plays no active role 
in our present simulations. On the other hand, the calcu- 
lations described in Sect. 14.71 demonstrate that there is a 
non-negligible amount of CO present in the lower chromo- 
sphere. Hence, a full treatment of CO as a cooling agent 
might even amplify the thermal bifurcation of the chro- 
mosphere ( see. e.g..lAvres|l98ltlAnderson fc Athavlll989t 
ISteffen fc Muchmorelll988h . 

Although the mean chromospheric temperature of our 
simulation lies considerably below the radiative equilib- 
rium temperature, we find a net radiative cooling of the 
chromospheric layers: (V • F ra( i}t is positive here. This ap- 
parent contradiction can be explained by the presence of 
sufficiently strong temperature fluctuations and the highly 
non-linear temperature dependence of the radiative heat- 
ing/cooling rates. We do not claim, however, that this 
situation is actually realized in the solar chromosphere. 

Our simulations indicate that the wave generation is 
mainly controlled by the large-scale dynamical evolution 
of the granulation pattern (see Sect. 14.2)) . This is at vari- 
ance with the c l assical picture of the Lighthill-Stein theory 
l|Lighthilil952HSteirll9"67lll968|) where small-scale turbu- 
lent eddies make the main contribution to the acoustic en- 
ergy flux. Applying th e Lighthill-Stein theory to the Sun, 
iMusielak etaD l)l994h find that the acoustic flux spectrum 
shows a maximum near v > 15 mHz, and hence is domi- 
nated by "short period waves" . 

The presented simulation can marginally resolve tur- 
bulent eddies in the convection zone with wavenumbers 
up to fc max k 27r/(5 Air). According to the classical the- 
ory, eddies with wavenumber k mostly contribute to the 
acoustic wave spectrum at frequency uj — ku^, where 
Uk is the turbulent velocity of eddies with wavenum- 
ber k. Since Si 1 km/s, the simulation cannot de- 
scribe the turbulence spectrum beyond w max ^ 30 mHz, 
f max < 5 mHz. We conclude that our present model can- 
not resolve the small-scale turbulence which is responsi- 
ble for the sound generation in the Lighthill-Stein theory. 
The acoustic flux resulting from our simulation decreases 
monotonically with frequency, and so has little in com- 
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mon with the spectrum predicted by the Lighthill- Stein 
spectrum. Hence, it appears doubtful whether the classi- 
cal theory, based on the assumption of isothermal, homo- 
geneous, and isotropic turbulence, captures the essential 
physics of the violent, highly anisotropic layers at the top 
of real stellar convection zones. We argue that our numeri- 
cal simulation correctly represents the basic mode of wave 
generation, even at the present spatial resolution. 

6. Conclusions 

Based on a detailed 3-D simulation of the solar granulation 
and the overlying atmosphere, we have studied the gen- 
eration of waves by the time-dependent convective flow, 
and the wave propagation and dissipation in the higher 
layers. The most important improvements compared to 
previous numerical simulations are (i) self-consistent dy- 
namics without a need for a driving piston like done by CS 
and (ii) a high spatial resolution which is obviously nec- 
essary for modelling the small-scale structure of the solar 
chromosphere. On the other hand, the LTE treatment of 
the thermodynamics and the radiative transfer is certainly 
unrealistic in the chromospheric layers. We have presented 
evidence that some of the basic features seen in our model 
are nevertheless representative of the (non-magnetic) in- 
ternetwork regions of the solar chromosphere. 

The main result of the present investigation is the 
discovery of a complex network of hot filaments pervad- 
ing the otherwise cool chromospheric layers. Caused by 
interaction of standing and propagating hydrodynamic 
waves of large amplitude, the model chromosphere is a 
highly dynamical, spatially and temporally intermittent 
phenomenon. Its temperature structure is characterised 
by a thermal bifurcation: hot and cool regions co-exist 
side by side. Temperatures in the hot filaments are high 
enough to produce chromospheric emission lines, and the 
cool "bubbles" are cold enough to form molecular fea- 
tures. Thus, the chromosphere is hot and cold at the same 
time. This picture of the 3-D structure of the solar chro- 
mosphere has the potential to explain the apparently con- 
tradictory observational diagnostics which cannot be un- 
derstood in the framework of one-dimensional theoretical 
or semi-empirical models. 

The presence of strong spatial and temporal tempera- 
ture fluctuations has a remarkable consequence: the tem- 
perature minimum and the outward directed temperature 
rise inferred from semi-empirical models might be artifacts 
in the sense that they do not necessarily imply an increase 
of the average gas temperature with height. Our model 
suggests that the radiative emission can be sustained by 
the hot propagating shock waves even though the main 
fraction of the chromospheric layers is cool and the mean 
gas temperature profile shows an almo st monotonic de- 
crease - a conclusi on already reached bv lCarlsson fc SteirJ 
l)l994 11995L Il997j) on the basis of one-dimensional hydro- 
dynamical simulations. 

We conclude that improved 3-D radiation hydrody- 
namic simulations of the kind presented in this work 



are likely to lead the way towards a consistent physical 
model of the thermal structure and dynamics of the non- 
magnetic solar chromosphere which eventually can explain 
the various observational diagnostics. 
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